

# ===========================================================
# =                Libraries                                =
# ===========================================================
library(foreign)
library(lme4)
library(rjags)

# ===========================================================
# =                Data                                     =
# ===========================================================

dat <- read.dta("data.dta")


# ===========================================================
# =         Estimation using MCMC via JAGS                  =
# ===========================================================

ones <- rep(1, dim(dat)[[1]])
X <- with(dat, cbind(ones,age,inclow,inchi,olead,lright,male,tenurez,tradez,gdpz,inflz))
K <- ncol(X)
cntry <- rep(NA, dim(dat)[[1]])
j <- 1
for (i in sort(unique(dat$cntry))){
	cntry[dat$cntry==i] <- j
	j <- j+1
}
J <- length(unique(dat$cntry))
N <- dim(dat)[[1]]

data <- list("y"= dat$support, "X"=X, "cntry"=cntry, "N"=N, "J"=J, "K"=K)

load.module("glm")

model <- "model {
for (i in 1:N){
	y[i] ~ dnorm(y.hat[i], tau.y)
	y.hat[i] <- alpha[cntry[i]] + inprod(beta[],X[i,])
}

tau.y ~ dgamma(0.001, 0.001)
var.y <- 1/tau.y

for (k in 1:K){
	beta[k] ~ dnorm(0, 0.001)
}

for (j in 1:J){
	alpha[j] ~ dnorm(0, tau.a)
}
# mu.a ~ dnorm(0, 0.001)
tau.a ~ dgamma(0.001, 0.001)
var.a <- 1/tau.a
}"

write(model, "model2.bug")


# Initialize two chains from different starting values
# I explicitely set the RNG seed here to allow for complete reproducibility
setinits <- function(x) {
	list(
		list(beta=rnorm(11), alpha=rnorm(14),
    ".RNG.name"="base::Wichmann-Hill", ".RNG.seed"=1234),
		list(beta=rep(0, 11), alpha=rep(0, 14), 
    ".RNG.name"="base::Wichmann-Hill", ".RNG.seed"=4321))
}


m2 <- jags.model(file="model2.bug", data=data, n.chains=2, 
	inits=setinits())

mon <- c("beta", "var.y", "var.a")
update(m2, 10000)
m2.out <- coda.samples(m2, variable.names=mon, n.iter=20000, thin=5)

colnames(m2.out[[1]]) <- colnames(m2.out[[2]])  <- 
 	c("intcpt","age","inclow","inchi","olead","lright","male","tenurez",
		"tradez","gdpz","inflz","var.a","var.y")

geweke.diag(m2.out, frac1=.1, frac2=.5)
heidel.diag(m2.out)
gelman.diag(m2.out)
autocorr.plot(m2.out)
summary(m2.out)



